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Abstract. 

If deleterious mutations near a fitness maximum in a spatially distributed 
population are sufficiently frequent or detrimental, the population can undergo a 
fitness collapse, similarly to the Muller’s ratchet effect in well-mixed populations. 
Recent studies of one-dimensional habitats (e.g., the frontier of a two-dimensional 
range expansion) have shown that the onset of the fitness collapse is described by 
a directed percolation phase transition with its associated critical exponents. We 
consider population fitness collapse in three-dimensional range expansions with both 
inflating and fixed-size frontiers (applicable to, e.g., expanding and treadmilling 
spherical tumors, respectively). We find that the onset of fitness collapse in these two 
cases obeys different scaling laws, and that competition between species at the frontier 
leads to a deviation from directed percolation scaling. As in two-dimensional range 
expansions, inflating frontiers modify the critical behavior by causally disconnecting 
well-separated portions of the population. 
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1. Introduction 

Typical laboratory evolution experiments, especially in microbial populations [1], are 
conducted in a well-mixed environment. Theoretical treatments of evolution also often 
focus primarily on the well-mixed case [2]. However, in nature, populations will usually 
have some spatial variation. The spatial distribution of genotypes is particularly 
important in range expansions, in which populations invade new territory [3], For 
example, cancer cells may form solid tumors which invade healthy tissue [4], microbial 
colonies may spread over a surface such as a Petri dish and animals like butterflies 
and birds settle new territory over time [7j. A key impact of spatial variations on 
evolution is an enhanced noise (called genetic drift) due to number fluctuations in local 
regions with small effective population sizes at the frontier. This genetic drift can lead 
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to fluctuation-driven phase transitions, characterized by the extinction of a particular 
strain within a population. This paper focuses on fitness loss in populations due to such 
phase transitions in three-dimensional range expansions. 

Previous work has focused on evolutionary dynamics in populations with one¬ 
dimensional frontiers [SUaEUID]. The spatial dynamics and geometry of the population 
strongly influences its evolutionary dynamics, creating dramatic changes relative to a 
well-mixed population. For example, Krug and Otwinowski [8j have found that spatial 
fluctuations can lead to a fitness collapse that is reminiscent of Muller’s ratchet in well- 
mixed populations mm, in which the average fitness of the population declines as 
its fittest members go extinct via the combined effects of genetic drift and deleterious 
mutations. Effectively one-dimensional habitats may be realized, for example, at the 
frontier of a two-dimensional range expansion, such as a microbial colony on a Petri dish 
mm- In such expansions, growth is often confined to a thin layer at the population 
frontier. 

In this paper, we study the evolutionary dynamics at the edge of three-dimensional 
range expansions such as tumors, spherical microbial colonies grown in soft gels on 
liquid media, etc. We will compare expansions of populations experiencing deleterious 
mutations at non-inflating fronts, as in figure [[[a), with those living on curved, inflating 
frontiers, as in figure [[[b). In both cases, we focus on expansions with a single-ccll- 
wide frontier of active growth. This limit maximizes the effects of genetic drift and is a 
realistic model of, e.g., tumor populations which may have thin actively growing regions 
ra. To make contact with non-equilibrium statistical physics literature mm , it is 
convenient to refer to such expansions as d = 2 + 1-dimensional expansions (or d = 1 + 1 
for two-dimensional expansions with thin frontiers), where the 2 is the effective frontier 
dimension, and the +1 refers to the “time-like” direction in which the expansion spreads, 
as shown in figure [Q 

Deleterious mutations often accumulate in genetically unstable populations, such 
as tumors with a high mutational load [16] or in viral populations, with their 
imperfect proofreading machinery in. Understanding the nature of this accumulation 
is important, for example, for designing effective cancer therapies [16]. Since the target 
space for deleterious mutations is typically much larger than for beneficial ones, there 
may be time scales over which populations almost exclusively experience irreversible, 
deleterious mutations. Hence, we will focus on the accumulation of such effectively one¬ 
way mutations and the onset of fitness collapse that can result. As discussed below, the 
collapse can be associated with a phase transition, slightly smeared by finite size effects. 
Much like at a second-order equilibrium phase transformation, we expect that universal 
behavior appears near this transition, with results insensitive to the microscopic details 
of our models [HJ. Indeed, many diverse systems can share a single universality class and 
identical scaling characteristics near the transition. For example, Otwinowski and Krug 
Pj showed that the onset of fitness collapse at the one-dimensional frontier of a flat (i.e., 
non-inflating) two-dimensional range expansion is governed by the directed percolation 
(DP) universality class, which describes a wide variety of systems, ranging from growing 
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Figure 1. (a) A small simulated three-dimensional range expansion (with periodic 

boundary conditions perpendicular to the time-like growth direction) and an initial 
population of all red cells, representing a maximally fit lineage. As the cells divide and 
populate the subsequent generations (the two-dimensional sheets of cells in a triangular 
lattice along the vertical direction) they acquire deleterious mutations (each decreasing 
the growth rate by 10%) with a 20% probability per division. The cells gradually turn 
more blue as they acquire more mutations, (b) A simulated spherical range expansion 
with the same deleterious mutations, cut in half to reveal the initial red population 
(about six cells in diameter) at the center. The dashed black lines indicate the edges of 
the two hemispheres. The generations are concentric spherical shells about a single cell 
wide. The population frontier inflates in such an expansion: The initial ball of red cells 
grows to a cluster with about a 20 cell diameter. Section [2] discusses the simulations 
in more detail. 


interfaces to turbulent liquid crystals pa E]. We describe here the analogous onset in 
three-dimensional range expansions and find important differences, including a different 
phase diagram shape and deviations from conventional DP scaling. We also consider 
the effects of an inflating frontier [see figure Pb)], which leads to a cutoff in the critical 
scaling by mitigating some of the effects of genetic drift. 

We model range expansions in which we track M different types of cells, where M 
may be finite or infinite. Each type will represent a cell with m accumulated deleterious 
mutations, where m = 0,1 ,... ,M — 1. We label the species with different colors, 
as shown in figure HI The red cells are the fittest cells with m = 0. The cells turn 
more blue as they acquire more mutations. For flat frontiers, as in figure Pa), we will 
study how increasing M influences the onset of fitness collapse. To study the effects of 
inflating frontiers, as in figure Pb), we will focus on the M = 2 case. An experimental 
realization of our models is within reach: Current advances in fluorescence detection have 
greatly expanded our ability to track multiple lineages of cells in spatially distributed 
populations using confocal fluorescence microscopy. In fact, tracking M fa 100 lineages 
using as many different fluorescent reporters may be feasible [18 ] fl9 l [20]. As such tools 
become more widely used, it is of interest to map out the theoretical space of possibilities 
for the evolutionary dynamics. 

This paper is organized as follows: In section [2] we introduce our range expansion 
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lattice model and a coarse-grained description using stochastic partial differential 
equations. Section [3] presents a mean-field analysis of the coarse-grained model which 
recapitulates some known results for well-mixed populations (such as the quasispecies 
model of Eigen and collaborators (2Tj) and illustrates various possibilities for the critical 
behavior. We present results for the onset of a fitness collapse in three-dimensional range 
expansions for large values of M with flat fronts [as in figure [D(a)] and compare to two- 
dimensional expansions in section [4] We focus on changes due to inflating frontiers for 
M = 2 in section [5] and conclude in section [6] 

2. The Model 

The expansions in figure |Tj were simulated using a variation of the models presented in 
PI [22]. In these expansions, deleterious mutations accumulate at a constant rate // per 
cell division. Each mutation has a multiplicative fitness cost, and the growth rate T m 
of a cell with m acquired mutations is thus 

r m = (i - 8 ) m , (l) 

where s (0 < s < 1) is the strength of the deleterious effect, and the “fitness class” m 
ranges from 0 to M — 1, so that M is the total number of cell types. Note that T m is 
always positive, and that the growth rate when no deleterious mutations are present has 
been rescaled to unity. In the small s limit, our model amounts to a spatial generalization 
of a well-mixed population genetics problem studied in [23]. To simulate planar fronts in 
a three-dimensional expansion, we stack two-dimensional “sheets” of cells in a triangular 
lattice (with periodic boundary conditions) to form a three-dimensional hexagonal close 
packed (hep) lattice. Each two-dimensional sheet of cells represents a single generation 
[see figure [T][ a)]. Each approximately square sheet will have N cells arrayed along the 
two dimensions of the sheet. At the frontier, the structure of the hep lattice ensures that 
each empty site is adjacent to three frontier cells, all of which compete to divide into 
that site, as shown in figure [2]( a). These three frontier cells divide with a probability 
proportional to their growth rate. In particular, a cell of type m is placed in the 
empty site with a normalized probability p m oc N m T m , where N m is the number of 
competing frontier cells with class m. Upon invoking a normalization condition, we 
have p m = N m T m / After a new daughter cell is placed in the empty spot, 

it can acquire an additional deleterious mutation moving it down a fitness ladder with 
some probability p as shown in figure [2](b). 

All empty sites in a two-dimensional sheet are filled before moving on to the next 
generation, ensuring flat population fronts for the three-dimensional range expansions. 
As in [3j, a single triangular lattice sheet may be used to simulate a two-dimensional 
expansion. In this case, rows of the triangular lattice represent the successive population 
frontiers and two cells compete for each empty site. Such uniform fronts are reasonable 
approximations when the population has an effective surface tension, such as at the 
periphery of a yeast colony [23], and when differential growth rates between species at 


Fitness collapse in spatial population genetics 


5 


the frontier are small (i.e., s -C 1 in our model). However, when present, undulations of 
the population frontier can strongly modify the evolutionary dynamics. For example, in 
two-dimensional expansions, deviations from directed percolation scaling were found at 
the onset of fitness collapse f25j. Although we do not study them here, undulations could 
cause deviations in the critical exponents associated with three-dimensional expansions, 
as well (see, e.g., TO- 

= 3 

n r 

_ iy m*- m 

~ ZfXeTe 

cells at frontier 

(b) 

m = 0 1 2 ... M -1 

Figure 2. (a) A single update step in our lattice model of planar range expansions. 

Cells at the frontier compete to divide into the (uncolored) adjacent empty site. 
The site gets filled with a cell of fitness class m with probability p m , determined 
by T m = (1 — s) m and N m , the number of competing cells with class to. In this case, 
three total cells are competing, so Ng = 3. (b) After this competition, the daughter 
cell can acquire an additional one-way deleterious mutation with probability /z, moving 
it from class to. to to + 1. The most deleterious class is to = M — 1. 

Special techniques are required for circular or spherical expansions to avoid lattice 
artifacts p]. For spherical range expansions, we pick an empty site closest to the center 
of the growing population that neighbors at least one cell along the population frontier. 
This update rule creates a uniform front and mimics a surface tension in the spherical 
cell cluster m The site positions for a disordered lattice are generated in advance 
using the Bennett algorithm, originally used to model disordered metallic glasses [27] 
(see 1221 for more details). Initial conditions in such expansions are specified by assigning 
particular colors to all cells a distance of Rq or less from the population center. In our 
simulations, Ro is measured in cell diameters a. After the expansion grows out to some 
maximum radius, we break up the resulting ball of cells into concentric spherical shells 
that are a single cell thick. These shells then define a temporal succession of population 
frontiers. Thus, in both types of expansions illustrated in figure [U the update rules are 
designed so that the population expands outward by one cell diameter a per generation. 
We measure time in generations and choose length units that fix the population front 
speed to be v = a/r g = 1, where r g = 1 is a generation time. 

We will be interested in the deleterious one-way mutation regime (s, p > 0) where 
we find a fitness collapse transition. We will not consider beneficial mutations (s < 0), 
which have their own interesting phenomenology [8]. Our mutations are irreversible, so 
that when the fittest class rn — 0 dies out, it cannot be regenerated. The system is then 
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in an absorbing subspace of states from which it cannot escape. Note that with this 
model we can have irreversible extinctions at all class levels m as long as all cells in the 
fitter classes n < m have died out as well. When M = 2, our lattice model for range 
expansions with planar fronts reduces to a d = 2 + 1-dimensional version of the Domany- 
Kinzel stochastic cellular automaton [28j. In this case, there is a single absorbing state 
in which the entire population is in the less fit m = 1 class. Transitions into a single 
absorbing state, governed by the directed percolation (DP) universality class, are the 
subject of a wealth of literature in nonequilibrium statistical dynamics m- Hence, as 
in the original d = 1 + 1 model (see 0 ). we expect that when p,s > 0, there will be 
a line of DP phase transitions at which we find the onset of fitness collapse PEUE]. 
We will use the M = 2 case to guide our thinking about the more complicated and 
biologically relevant limit M > 1. 

To find an analytically tractable description of the evolutionary dynamics, it is 
useful to move to a coarse-grained model. Let f m = f m (x.,t) be the fraction of m-class 
cells in a small, approximately planar patch surrounding a point x on the population 
frontier at time t [9]. The cells in each local patch are treated as if they are in a 
well-mixed population. The cells are allowed to hop to an adjacent patch via a local 
cell rearrangement (occurring, for example, via outward pushing due to cell division or 
slight cell motility). The evolutionary dynamics in each well-mixed patch is modelled 
by a stochastic differential equation equation, developed by Good and Desai [23], that 
assumes the growth rates in equation (P) in the small s limit, a constant population 
density, irreversible mutations between classes m and m +1 with rate p, and the genetic 
drift associated with the finite population size in each patch. Then, for a non-inflating, 
flat population frontier, the fractions f m of m-class cells obey the spatial generalization 
of this equation: 

d t fm = DV 2 /m + s((m) - m)f m - nf m ( 1 - 8 m>M - i) + h/m-i 

M—l 

+ - /m)y / 2r- 1 A/ £ ^(x,t), (2) 

t=o 

where is a Kronecker delta function = 1 if i = j and 0 otherwise), (m) = 
S^c) 1 * s the local average fitness class, and r g is the generation time. The strength 
of the genetic drift, A = A -1 , is set by the population size N of each local patch. For 
each class, the Ito noise terms &(x, t) [29] have vanishing mean and variances given by 
(&( x , t)£ m (x', t')) — Si m S(x — x')S(t — t'). The diffusion constant D ~ a 2 /r g describes 
random cell rearrangement at the frontier, which in our lattice simulations is expected to 
move cells by approximately a cell diameter a per generation. We impose the boundary 
conditions /_i = 0 and ft — 1- Our model is essentially a stepping-stone model 

of the range expansion M E], and can be derived more rigorously from a modified, 
continuous-time version of the lattice model [9]. Note that our lattice simulations are 
in the limit where each patch population size is N = 1, which maximizes the effects of 
genetic drift. 
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We will examine equation ((2]) near the fitness collapse transition, where the fractions 
f m at the fit edge of the fitness distribution (i.e., the first few values of m) will be small 
and close to their absorbing state values f m = 0. When the fractions f m are small, 
equation (J2]) reduces to a set of (nonlinear) coupled DP Langevin equations for the 
first few values of m, i.e., directed percolation with “colors,” studied extensively in 
[3T1 [32], [33]. To see this, we set Jm-i = 1 — Y1?c=q 2 ft an d keep on ly the leading order 
noise terms in equation ((2]) . Then, the Langevin equations for 0 < rn <C M — 1 reduce 
to: 


d t fm ~ D7J 2 fm — 


M—2 

T m E Tmtfe 

1 = 0 


fm + Ffm-l + ^Tg 3) 


where r m = //— s(M— 1— m) and T m i = s(M—1—£). One can show that the higher order 
noise terms we dropped are irrelevant in the renormalization group sense. Equation (j3J) 
is now a hierarchy of DP Langevin equations, coupled to each other both linearly, via 
the mutation term nf m - 1 , and bilinearly, via the off-diagonal matrix elements 
(with m 4 i) |3T, [321 S3]- When M = 2, equation (j3J) reduces to the DP Langevin 
equation [14] for the m — 0 class, studied in the context of spatial population genetics 
in [9j. Hence, for the special case M — 2, we expect to find DP critical behavior in 
our model. We will confirm this result in the next section. It is possible to generalize 
equation (J3J) to include the effects of inflating frontiers [221 EJ. Although we study 
inflation via simulations using a disordered lattice in section [5] an analytic treatment of 
the inflationary generalization of equation (J3J) is beyond the scope of this paper. 

If we allow r m and T m i in equation (j3j) to be arbitrary (instead of given by the 
expressions just below equation (J3])), this set of equations exhibits a rich phase structure. 
In particular, the parameters {r m } can be tuned so that any, some, or all of the m-species 
go extinct. The multicritical point at which all of the species go extinct simultaneously 
(for M —>■ oo) is particularly interesting as it describes the critical dynamics of some 
interface growth models at a roughening transition [31] [8]. The held / m (x, t) for the 
interface models represents the fraction of interface plateaus of height m at some point 
x and at time t. A simple Langevin equation exhibiting such a multicritical point 
was proposed in ra- The equation has the same form as equation (131) . but with 
constant coefficients r m = r and a constant, diagonal matrix T m £ = TS m £. This equation 
describing the interface growth models falls into the unidirectionally coupled directed 
percolation (UCDP) class near the multicritical point [33- 34]. We shall see that this 
UCDP universality class will be relevant for our range expansion evolutionary dynamics, 
as well. The presence of the bilinear couplings T m i with m 4 £, however, can lead to 
deviations from unidirectionally coupled directed percolation scaling [35j. We will now 
describe a mean-held analysis of equation (J3J) to better understand the phase structure 
of our model and to point out various possibilities for the critical behavior. 
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If we remove the noise (i.e. genetic drift) and diffusion terms in equation ([3]) by setting 
D — A = 0, we recover the mean-field limit. This limit also corresponds to the dynamics 
of an infinite (N —> oo), well-mixed population. In this case, our fitness class fractions 
/ m (x, t ) are independent of x and evolve deterministically. Hence, we cannot have noise- 
driven extinction events. However, there is still a phase transition for any finite total 
number of species M < oo. Since this transformation is not noise-driven, it is not 
a Muller’s ratchet and is more properly called an “error threshold” transition [36], a 
term originating from the theory of quasi-species At and above this threshold, the 
mutation rate is high enough to induce a fitness distribution collapse in which the entire 
population eventually collapses into the least fit m — M — 1 class (largest value of m), 
as illustrated in figure [3] 



Figure 3. An error threshold transition for a well-mixed population. When 
g/s < M — 1, we are in an “active” phase where the population has a fixed fitness 
distribution at long times, determined by mutation-selection balance. The distribution 
is given approximately by equation (j8]), i.e. a Poission distribution with mean g/s 
(lower inset). When g/s > M — 1, the fitness distribution collapses and the entire 
population is eventually in the least fit class for long times. 


We can solve for the mean-field stationary fractions —>• oo) using 

the method of generating functions EH- The asymptotic long time fractions f^°° obey 
the steady-state, mean field version of equation (J3]): 




M—2 


Tmlf 1 


t—>oo 


1=0 


nt^fOO m 

J m '“M 


(4) 


with r rn — p — s(M — 1 — m) and T m (, = s(M — 1 —£). The generating function packages 
information about the different fitness classes into a polynomial function of x, defined 
by 


M—l 


G(x) = V /P°V. 
e=o 


( 5 ) 
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Equation flj]) simplifies to a differential equation for G(x) when we multiply both sides 
by x £ and sum over all possible £: 


dG 

dx 


-G(x) + 
s 


a 

s 


1 

X 


1 


rt—^oo M—1 
JM- 1 x j 


( 6 ) 


with the boundary condition G(l) = 1. The last term in equation (]6jh proportional to 
the fraction of cells in the last allele class, f/jFf , can be reexpressed in terms of G(x) 
using the relation 


/, 


£->-oo 


d r ‘ 


m\ dx m 


G(x) 


x=0 


(7) 


Upon solving equation (j6]) for G(x), we use equation (J7J) to find the steady-state fitness 
distribution. Equation (J6]) is most readily solved when the fraction fjyFf is small, or if 
M —> oo. Then, we can ignore the term proportional to in equation (El) and find 

that G(x) = which is the generating function for a Poisson distribution: 


/y oo = T(-)’V'‘ / ‘ (/«“ « 1 or M -+ oo). (8) 

Equation ((8]) will be a good approximation to the steady-state distribution as long as 
p/s <C M — 1 in a well-mixed population, i.e. far below the dashed transition line 
in figure |3l For spatial range expansions, however, the steady-state distribution will 
be quite different. To highlight this difference, we will calculate f (the fraction of 
the population in the most fit class) for the range expansions in the next section and 
compare to the well-mixed population result from equation (18]) : /= e~^ s . 

When M < oo, equation (EJ) exhibits a line of phase transitions in the (s, /i)-plane 
at r 0 = p — s(M — 1) = 0, i.e. for p/s = M — 1. If r 0 > 0, the entire population is 
eventually forced into the lowest fitness class, so that f^°° = d m ,M- i • Conversely, if 
To < 0, we maintain a stationary fitness distribution spread out over many fitness classes. 
The phase diagram and steady-state solutions are illustrated in figure El A good order 
parameter for the phase transition (which will also be useful when we analyze range 
expansion simulations) is the steady-state fraction of cells in the fittest class, / 

If f^°°{p, s) = 0, we are then in the inactive phase where the cells in the m — 0 
class cannot survive at long times due to the deleterious mutations. Conversely, if 
ft*°{p, s) > 0, we are in the active phase and the most fit m — 0 cells constitute a 
non-zero fraction of the population, represented by the left edge of a Poisson distribution 
(see lower right inset of figure E])- This phase transition in the infinite population limit is 
not unique to our model. For example, an analogous fitness distribution collapse occurs 
in an infinite-population model of the evolution of DNA binding sequences [ 38] . 

When M — y oo, the order parameter is non-zero for any p, s > 0, since f/^ 00 = 
e~ ,J ‘/ s > 0, and the phase transition disappears! However, this behavior is an artifact 
of the limit of infinite population size. If we had some finite population size N and 
introduced the noise term back into equation (13]), then a transition can arise if the 
number of individuals in the rn — 0 class, no ~ Ne j s driven to zero via number 
fluctuations, leading to an irreversible loss of individuals in the highest rn — 0 fitness 
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class, and a shift in the fitness distribution. A succession of such shifts is known as 
Muller’s ratchet in a well-mixed population [12], The differences between Muller’s 
ratchet and the error threshold in a well-mixed population has been discussed further in 
population genetics literature [39] . The lack of a transition for M-h» arises because 
Muller’s ratchet does not operate in a noise-free infinite, well-mixed population. The 
finite-M error threshold transition is the relevant mean-held description of the fitness 
collapse in the spatial range expansions of interest here: it accurately describes range 
expansions above the upper-critical dimension d = 3 + 1. However, in biologically 
relevant dimensions, such as the d — 2 + 1 and d — 1 + 1 cases discussed in the 
introduction, the diffusion and noise terms in equation ([3]) significantly modify the 
mean held result. 

Consider the mean-held transition in more detail: When we approach the critical 
point from the active phase, r 0 —> CD, the expression r m = p — s(M — 1 — m) simplifies: 
r m —> sm > 0 for all m > 0. Thus, r m is positive, and all of the other classes m > 0 
are already in their inactive state. Consequently, the cells with classes m > 0 have 
no interesting critical dynamics of their own and inherit the critical scaling of the 
m — 0 class via the mutation term pf 0 . So, as r 0 —> CD, the steady-state fractions 
of all the species vanish according to ~ |ro|^ MF (0 < m < M — 1), where the 

critical exponents are (5^ = 1. In other words, the critical exponents assume the m- 
independent value given by the mean-held DP exponent governing the critical behavior 
of the m = 0 class. Thus, the densities of all fitness classes vanish in the same way 
within mean-held theory. 

In d = 1 + 1 and d = 2+1-dimensional range expansions (and even for d = d'+l with 
d' > 3), the noise and diffusion terms in equation (|3]) will renormalize the coefficients r m 
and T m £. These renormalizations lead us to consider general couplings r m and T m £, which 
leads to a rich variety of possible critical behaviors. For example, consider the M —)■ oo 
case. Otwinowski and Krug argued that the transition to a fitness collapse in a two- 
dimensional range expansion (d = 1 + 1) with M —> oo will fall into the same universality 
class as the interface growth models which exhibit multicritical, unidirectionally coupled 
directed percolation (UCDP) behavior [8[|3T]. The UCDP scaling regime is achieved in 
equation (J3]) when we hx T m e = T5 m i and let r rn — r —> 0“ for all m > 0 [33.- [34] 31 . 
In this case, the exponents (3^ for each class m are quite different! In particular, by 
analyzing equation (J4]) with these new coefficients, one finds that ~ |r|^MF as 

r —> 0+ where (3 ^p = l/2 m [31]. Note that in this case, only the m = 0 class fraction 
has mean-held DP scaling. The other fractions have exponents that rapidly 
decay with m. This is the mean-held version of the unidirectionally coupled directed 
percolation universality class. If we include the off-diagonal terms in T m we may find 
universality classes different from UCDP. These off-diagonal terms represent competition 
between different classes and become more relevant in higher dimensions [35]. The lattice 
simulation results in the next section strongly suggest that the competition terms modify 
scaling in three-dimensional range expansions when M —> oo. 

Another way to describe the htness collapse transition is to study how the fractions 
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fm it) decay with time at the transition. Ideas from directed percolation theory lead 
us to expect that these fractions vanish according to the power laws f m (t) ~ t~ s(rn) 
where S ^ /o^' n> [13]. Here, the exponents are time-like correlation length 

exponents for each class m (see next section for more details). In the mean-held 
approximation, = 1 and = 1 f° r all m dr equation ([3]). At 

the mean-held UCDP point, however, = l/2 m , which decreases rapidly with m. 
A priori, it is not clear which (if any) of these mean-held behaviors describes the 
evolutionary dynamics of range expansions, since the diffusion and noise terms are 
expected to modify the critical behavior. However, the mean held results help map out 
the space of possibilities. We will use simulations in the next section to calculate the 
exponents m m ) for two- and three- dimensional range expansions for the hrst few fitness 
classes m. 

In addition to changing the scaling exponents, the noise and diffusion terms in 
equation (|3]) shift the critical line in the (s, / u)-plane relative to its mean-held counterpart 
r 0 = /a — s(M — 1) = 0. We expect that the local genetic drift at the population frontier 
will induce a transition even for M —> oo. In the Krug and Otwinowski model with 
M — y oo, for example, the transition in two-dimensional range expansions occurs when 
p/s 2 ~ 1 [8]. We shall see in the next section that there are important differences 
for three-dimensional range expansions. Although there is a fitness collapse transition 
when M —> oo, as in two-dimensional range expansions, and we find indications of 
multi-critical scaling similar to the UCDP universality class, the exponents we find are 
significantly different from expected UCDP scaling for d = 2 + 1 [33]. Also, the phase 
diagram shape is qualitatively different from the two-dimensional case. 

4. Lattice Simulation Results 

As discussed above, a natural order parameter for the htness collapse transition is the 
infinite time limit of the fraction of cells in the fittest class, m = 0, averaged over many 
range expansions with identical initial conditions. This fraction, /g^°°(/i, s), is then 
also averaged over the entire population front. The initial condition is a population 
composed exclusively of cells in the most fit m = 0 class. We approximate s ) in 

our simulations by computing the fraction of m = 0 cells at some long time t at which 
the system has already settled into a steady-state. The corresponding phase diagrams 
in the (s, /i)-parameter space for d = 1 +1 and d = 2 + 1 are shown in figure 0J Near the 
fitness collapse transition, the system takes longer and longer to settle into the steady- 
state due to critical slowing down. However, at the resolution of the phase diagram 
of figure [H the sampling time is long enough so that the corresponding corrections 
are negligible. The phase diagram we find stops changing with increasing M for both 
d — 1 + 1 and d = 2 + 1 when M > 10. Hence, to understand the M —> oo case with 
less computational effort we just look at M = 40 (as well as M = 2) for d = 2 + 1 in 
figure [4l For d = 1 + 1, we take M = 20 as a proxy for the limit M oo. 

The phase transition line shape for two-dimensional expansions has been calculated 
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Figure 4. The phase diagram illustrating the transition from a stationary fitness 
distribution (active phase) to a collapsing fitness distribution (inactive phase) for three 
and two dimensional range expansions. The heat map describes the case of a flat front 
with M = 40 for d = 2 + 1 dimensions, with M = 40 representing our approximation to 
the M —> oo limit. The color plot shows the fraction of cells fo in the most favorable 
fitness class at t = 2000 generations (averaged over 300 expansions) for an initial 
population frontier of N = 500 2 cells (a triangular lattice 500 cells long and 500 cells 
wide) all in class 0. We believe that t = 2000 closely approximates the steady-state 
value We find precise locations (black crosses) at which the transitions for 

M = 40 and M = 2 occur by checking for dynamical critical scaling. The transition 
line for M = 40 (solid line) is given approximately by g ~ 0.81s. The M = 2 transition 
has a prominent logarithmic correction and can be fit by g ~ 1.86s/ln(s/24). The 
transition in d = 1 + 1 dimensions (with M = 20 effectively approximating the M —> oo 
limit) has a very different shape, with g « 0.5s 2 , consistent with results of [8} |5]. 


previously M- It is given by p oc s 2 , where the constant of proportionality is model- 
dependent and will vary with M. The particular scaling /i oc s 2 , however, does not seem 
to change with M. In our lattice model, we find that p ~ 0.5s 2 for the large M case, as 
shown by the dashed line in figure [U To calculate the shape of the transition line for 
three-dimensional (i.e. d = 2 + 1) expansions, we start with the M = 2 case and then 
check with simulations that only non-universal parameters change with increasing M. 

The transition line shape for three-dimensional expansions with M — 2 may be 
calculated by analyzing the scaling near the special s = p = 0 point. This point is 
described by the well-known voter model ns at which the two strains compete neutrally 
with no mutations. A non-zero s along the line p = 0 introduces a bias to the voter 
model and controls a phase transition. Without mutations, our initial condition of all 
m = 0 cells is unchanged by the dynamics, and cannot be used to define an order 
parameter. Instead, we take as our p — 0 order parameter the survival probability of 
m = 0 cells starting from an initial population frontier of one m = 0 cell (at the origin, 
say) surrounded by m = 1 cells. For s < 0, the m — 0 cells engendered by this isolated 
mutant always die out (inactive phase), while for s > 0 there is a non-zero survival 
probability (active phase). The critical point itself is at s = 0. Note that the voter 
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model has its own universality class, different from directed percolation mm- 

A non-zero mutation rate p > 0 pushes the system out of the voter model class. The 
corresponding cross-over scaling description has been calculated previously using field- 
theoretic techniques [32]. The transition line shape, derived from cross-over scaling, has 
the form p = As/ | ln(s/5)|. Both the constant of proportionality, A, and the constant 
B inside the logarithm are non-universal. A similar transition line shape is found for 
the onset of mutualism in three-dimensional range expansions [26J. We find these non- 
universal parameters by fitting to simulation data, as shown in figure [H This cross-over 
scaling shape works well for M — 2, but the logarithmic correction is not as evident 
in the large M case (B is much larger than typical values of s). Hence, we fit the line 
shape to p — As in figure [4], with A as our single fitting parameter. This linear shape 
works extremely well (see M = 40, d — 2 + 1 transition line in figure [4]). 

Scaling near the voter model point also determines the behavior of /o _>0 °(yU, s) deep 
in the active phase, when /q^°°(/x,s) is close to 1. In a well-mixed population, this 
means p/s <C M — 1 and we are far away from the transition. As discussed in section [3l 
ps e~^/ s ps 1 — p/s in this case. In d — 1 + 1 dimensions, f/^ 00 may be calculated 
by mapping the boundaries of localized patches of less fit (m > 0) strains to random 
walks mmn- This mapping is essentially the same one employed to study clusters 
in the d = 1 + 1 voter model. However, the steady-state fraction now depends on 

the scaling combination /i/s 2 , instead of p/s. The simulation results for fo~*°°(p, s) are 
presented in figure [5](a) for d = 1 + 1 and d — 2 + 1. In d — 2 + 1, the mutant cluster 
shapes are quite complicated, as illustrated in figure [5](b). A simple mapping to random 
walkers is not possible in this case. However, we will now make a scaling argument for 
fo^°° that exploits exact voter-model results for d = 2 + 1. 

Consider that deep in the active phase, there will be finite-sized connected clusters 
of species with fitness classes m > 0, and in class m — 1 in particular. The rest of the 
population will be in the fittest m = 0 class. Provided clusters of cells with m — 1 
are so dilute that they do not collide, we can treat each cluster independently. The 
m = 0 state, in this case, acts as the absorbing state for each cluster. Note that this is 
backwards from the analysis of the fitness collapse transition, where we treat the m = 0 
state as the “active” state. The non-interacting rn — 1 clusters with some finite average 
width £j_ and length are shown in figure [5](b). If the m = 0 class has some selective 
advantage coefficient s > 0, the m = 1 clusters will have a relative selective disadvantage 
— s < 0. So, the m = 1 clusters will be governed by voter model scaling in the inactive 
state. They will contain some average number of cells (n) s , which will depend on s. 
Since the clusters do not collide, we can approximate f/^ 00 as 

i-ddi, (9) 

jVf r . 

where Af r . is the total number of cells at the population frontier. Also, since the cell 
populations with m > 1 will be very small, the value of f/^°°(p,s) should not vary 
much with M. 

We now estimate how the average number of cells ( n) s in an m = 1 cluster scales 
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Figure 5. (a) The average steady-state fraction of the most fit 0-class cells, fo(g, s), 
for many different values of g and s. For d = 2 + 1 and M = 40, the data set is the same 
as in figureQJa). For d = 1 + 1, we evaluated fo(g, s) at t = 10 4 generations (averaged 
over 2000 runs) for a frontier of size N = 2000 cells and M = 20 . We also included 
the <2 = 2+1, M = 2 case for comparison (N = 500 2 cells, t = 2000 generations, 128 
runs). Upon using the scaling variables in equation ED. the data points collapse onto 
a single curve near the region f 0 (g,s) « 1 that describes the active phase. For the 
d = 2 + 1 cases, we find A 3 d ~ 0.3, so ~ 40 for M = 40 and A 3 d ~ 0.46, so « 11 for 
M = 2. For d = 1 + 1, we find A_ 2 d ~ 1.05. (b) A snapshot of all the m = 1 clusters 
with M = 5 in a section of a d = 2 + 1 range expansion deep in the active phase 
(s = 0.08, g = 0.0008). These clusters an average width £j_ and length £||, illustrated 
schematically with the arrows. 


with s. Consider a cluster formed at the origin r = 0 at time t = 0. We may calculate 
the pair connectivity function Y(|r|,£; |s|), which is the probability of finding a cell in 
the cluster a distance |r| from the origin along the population frontier at time t [14], 
The scaling behavior of T(|r|, f; |s|) near the phase transition at s = 0 may be extracted 
from the general scaling relations derived in [02]. At d = 2 + 1, we are at the upper 
critical dimension of the voter model and there are logarithmic corrections. Using the 
general scaling relations in [02], we find that ( n) s satisfies 


(n) s ^ 


a T n 


d 2 r dt T(|r|, t; |s|) 
1 


a 2 T g X\ ln(A/A 0 )| 


d 2 r dt T 


\J A| ln(A/A 0 )| ’ ^1 ln(A/A 0 )| 


; A|s| 


( 10 ) 

( 11 ) 


where a is a cell diameter, r g a generation time, A an arbitrary scale factor, and Ao is 
some non-universal constant. The two-dimensional integral over r ranges over the entire 
population frontier. This integral is of order the frontier area, which is approximately 
Nf r a 2 . The |s(-dependence of ( n) s may now be extracted from equation ([IT]) using a 
judicious choice of A = 1/s and an appropriate rescaling of both r and t. After these 
manipulations, we find that in the inactive phase: 


ln(s/s 0 ) 


(n) s ~ A 3D N {r _ 


s 


( 12 ) 
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where A 3 d is a constant proportional to the residual integration in equation (fill) and 
s 0 = Aq 1 . Upon returning to our expression for f 0 in equation (J9]) , we find that 



(13) 


Upon combining this analysis with results for well-mixed and d = 1 + 1 dimensional 
systems, we expect the following small p behaviors in the active phase: 



well-mixed 


d = 1 + 1 


(14) 



where A 2 d, 3 d and s 0 are non-universal and depend on the details of the models. We check 
the validity of equation (jT4|) for d = 1 + 1 and d = 2 + 1 for small p with simulations in 


figure El For our model, we find A 2 d ~ 0.5, A 3 d ~ 0.3, and so ~ 40. These parameters 
yield good data collapses in the active phase, i.e., for values of x < 0.2 in figure [5] As we 
move away from the strongly active phase region, the approximations in equation (fT4l) 
start to fail due to cell cluster collisions. Eventually, we reach the fitness collapse regime 
which we now investigate in more detail. 


0.2 


0.15 



S Q 
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Figure 6. Effective exponents at the fitness collapse transition versus time t in 
a two-dimensional range expansion evolved for t generations with a one-dimensional 
frontier of N = 10 4 cells and an initial population entirely in the m = 0 class. We tune 
g to the fitness collapse transition for three different values of s: s = 0.05 (crosses) 
s = 0.1 (circles), and s = 0.2 (triangles). The lines connecting these points are to guide 
the eye. The black dashed lines show the expected UCDP exponents [34]. Previously 
calculated exponent values for m = 3 were not available. 


A signature of the phase transition, when p is tuned upward to reach fitness collapse, 
is the power-law decay of the fitness class fractions with time: 



( 15 ) 
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where 8^ = represents the family of critical exponents introduced in the 

previous section. A good technique for finding these exponents from simulations is to 
calculate the time-dependent effective exponent 8^ m \t)'. 


8^\t = U) 


ln[/m(*f)//m(*t-l)] 

lnfc/Vi] 


( 16 ) 


where the {ti} are the sampled times in our simulation, chosen so that U/U_i ~ 2. The 
exponent may be calculated by estimating the steady-state value of 8^ (t) at long times. 
The error may be estimated by looking at the fluctuations of 8^(t). This technique 
works well for the m — 0 class, but we must introduce a modification for m > 0. Just 
as in the interface growth models [31], the exponents converge faster if we look at the 
integrated fractions F m = Y1 T=o ft- Then, since f m -i -C f m near the fitness collapse 
transition, we expect that F m ~ t~ s(m> . We use these integrated fractions F m , instead of 
f m , to calculate 8^ via the effective exponent technique. We do this for both d = 1 + 1 
and d = 2 + 1. 

The effective exponents 8^ (t) for d — 1 +1 are shown in figure |6l calculated using 
the multiple color generalization of the Domany-Kinzel model discussed above (281 HI- 
We see that the exponents decrease with increasing m, as we would expect at 
the multicritical UCDP point discussed in the previous section. The effective exponent 
results are compared to exponents reported for UCDP in the literature (dashed lines 
in figure [6]) [3l] . Our results ~ 0.164(6), c^d ~ 0.087(12), c^d ~ 0.035(13)] 
are consistent with the UCDP exponents for m = 0,1,2. The m — 3 exponent does 

/'o'v 

not appear to have a reported value and we compute ^(t) ~ 0.012(8). Hence, the 
fitness collapse transition for two-dimensional range expansions appears to be governed 
by UCDP, the universality class of the interface models discussed in [31] [33], [34] . This is 
an important conclusion, because it means that the off-diagonal competition terms T m £ 
(m i) in equation ([3]) are evidently irrelevant in d — 1 + 1 dimensions. We shall see 
that the situation in three-dimensional expansions is quite different. 

In the inactive phase, there will be a local speed V associated with the collapsing 
fitness distribution. This speed is the analogue of the Muller’s ratchet “clicking” rate 
as the fitness diminishes in well-mixed populations. As discussed in [8], the exponents 
zzj m) govern the speed V. Previous results for d — 1 + 1 suggest that the exponents 
i/™' 1 in the UCDP class do not vary with m and are all equal to the DP value 

= u\\ pe 1.733847(6) [ll][3T]. So, in two-dimensional expansions, we expect that the 
fitness distribution speed is V ~ |rpi, where r = p/s 2 — n c /s 2 is the distance away from 
the critical point [with (s c ,p c ) a critical point on the line of DP transitions]. Finally, 
note that there is a slow, upward drift in the m > 0 effective exponents (^) as t 
increases in figure [6] This effect was noticed in other UCDP models [Mj. Hence, it is 
possible that the multicritical regime is only relevant at intermediate times and that 
all of the exponents eventually approach the DP value reached by the m — 0 effective 
exponent, as one might expect from mean field theory. Nevertheless, the multicritical 
behavior is clearly important for a wide range of times as the population evolves. 
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Figure 7. (a) The effective exponent ( t ) (defined as in equation (ITCl) ) describing 
the temporal decay of the fraction f 0 of 0-class cells in a three-dimensional range 
expansion with a frontier of N = 700 2 cells at the fitness collapse transition (averaged 
over at least 1000 runs). The range expansions have different values of M and an initial 
population of cells all in the most fit m = 0 class. The different sets of points (with 
connecting lines to guide the eye) at the same M correspond to different positions 
along the line of phase transitions shown in figure HJ For the M —> oo, the number 
of mutations a cell may acquire during a simulation run is unbounded, (b) Effective 
exponents for larger values of m for the M — > oo case. As in the d = 1 + 1 case 
(figure ini), the exponents decrease with increasing m , indicating a multicritical scaling 
regime. We tune /i to the fitness collapse transition for three values of s: s = 0.075 
(crosses) s = 0.1 (circles), and s = 0.15 (triangles). In both (a) and (b), the dashed 
line indicates the expected directed percolation (DP) exponent value « 0.451 ITfl . 

We now analyze the fitness collapse transition in three-dimensional range 
expansions. We’ve seen that the scaling combination of s and /j that determines the 
shape of the phase transition line is x = p/s (up to logarithmic corrections). Hence, 
there will be some critical value x c describing the transition line in the (s, /i) plane. The 
parameter r = x — x c will measure the distance away from the transition. For M = 2, we 
expect the 0-class fraction at r = 0 to decay with the DP exponent ~ 0.4505(10) [14]. 
The results for (t) for various M are shown in figure [7](a). As expected, the effective 
exponents f° r small M = 2,3,4 approach the DP value (dashed line). However, 

as M —» oo, we find a significant deviation from DP scaling. The deviation seems to 
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occur at around M > 5. We find a scaling exponent ^ ~ 0.53(2) for M —* oo that is 
higher than the expected DP value [circles in figure 0(a)]. This value does not appear 
to change along the transition line. Thus, our exponent deviates from what is observed 
in the interface growth models. Note that in the Langevin equation proposed for the 
interface model [[3lJ, the 0-class dynamics obey the DP equation [14]. Our equation, 
the m = 0 case in equation (J3J), is different because it includes off-diagonal terms in 
the matrix Tmi- These competition terms might be responsible for the deviation from 
DP scaling at large M. This deviation is also expected from mean-held theory, where 
adding the competition terms can increase the exponent <5^ by a factor of 2 [35] . 

We also check to see if there is multicritical behavior in this model when M —)■ oo. 
We do find that when M —^ oo, the higher order exponents 5^ decrease with increasing 
m, as shown in figure [7](b). The exponent values we find, however, are quite different 
from the expected UCDP values. We find & 0.53(2), ~ 0.38(2), & 0.24(2), 

and hgp ~ 0.14(3), compared to previously calculated values ~ 0.46(2), 

0.26(3), and ~ 0.13(3) for UCDP with d = 2 + 1 [34]. We also see an upward shift in 
the exponents over time, just as in the two-dimensional range expansion case in figure [6] 
Hence, the multicritical regime might be transient. Nevertheless, the m — 0 effective 
exponent, (^(f), appears to settle to a fixed value. So, there is clearly some interesting 
deviation from regular DP behavior as we increase the number of species M. Assuming 
the exponents for various m do not deviate from the directed percolation value, 
the speed of the fitness wave associated with the fitness collapse should be governed by 
d = 2 + 1-dimensional DP [14]: 

V ~ |rl^ 11 with v\\ ~ 1.2950(60). (17) 

Our simulation results are consistent with this scaling of the speed V (data not shown), 
but a thorough check of the ujj m ^ exponents is beyond the scope of this paper. More 
extensive simulations would be necessary to verify that our model falls into a class 
distinct from DP and is not exhibiting a long-lived transient. For biological applications, 
however, our analysis is relevant because the multicritical behavior can have measurable 
effects over thousands of generations. 

5. Effects of Inflation 

Spherical range expansions also exhibit a fitness collapse, as illustrated in figure [8] Deep 
in the active phase, the inflationary nature of the population frontier will be largely 
irrelevant because, as discussed in the previous section, the evolutionary dynamics will 
be governed by small clusters with m > 0, which are insensitive to the front inflation. 
Similarly, large enough mutation rates will eventually extinguish the 0-class individuals 
in the inactive phase. Inflating frontiers will most strongly influence the evolution near 
the onset of fitness collapse. As argued in an earlier study of two-dimensional, circular 
expansions [9], inflation will causally disconnect portions of the population and prevent 
correlations from propagating along the frontier. There is still an abrupt population 
collapse, but certain critical properties are modified, similar to a finite size effect. 
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active critical inactive 



Figure 8. Spherical range expansions also exhibit a transition, but certain features 
are smoothed out, similarly to a finite size effect. At different times t (measured in 
generations), we show the surfaces of growing spherical populations with M = 6 and 
initial radii of about ten cell diameters, Rq ~ 10a. The initial populations have all 
0-class (red) cells. The cells turn more blue as they acquire more mutations. We show 
the population in the active region [p = 0.05, s = 0.2), inactive region (p = s = 0.1), 
and near criticality (p = 0.1, s = 0.14). 


For simplicity, we study the modifications due to inflation for spherical expansions 
just for the M — 2 case. The comparison between inflating and non-inflating frontiers 
for M — 2 is particularly informative, because we know that the non-inflating expansions 
are characterized by a genuine DP transition. Higher values of M will also exhibit a 
transition, as illustrated in figure [BJ We expect the effects of inflation to be similar 
for increasing M (and the M —> oo case in particular), but a direct comparison 
would require a better characterization of the universality class of non-inflating, three- 
dimensional range expansions. As discussed in the previous section, the M —* oo case 
appears to violate DP scaling, but a thorough investigation of this potentially new 
universality class is beyond the scope of this paper. 

An important confounding factor is the presence of a finite system size. Note that 
in the spherical expansion illustrated in figure [TJb) and figure O the genetic patterns 
evolve on a finite spherical surface. Simplifications arise for flat range expansions like 
those in figure [0(a), where an arbitrarily large frontier area mitigates finite size effects. 
To disentangle the effects of the finite surface area and the effects of inflation, we will 
compare inflating spherical expansions, as in figure EJ to populations at the surface of 
treadmilling spheres of some fixed radius Rq. The cells at the frontier of these fixed 
radius spheres will divide and displace each other over time, thus creating a dynamical 
“trcadmilling” effect. As discussed in the introduction, such a treadmilling sphere can 
model an avascular tumor which turns over cells at its surface but does not expand due 
to apoptosis at its center and pressure from the surrounding tissue [13]. The treadmilling 
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sphere simulations are performed by growing an initial spherical population of radius 
R 0 out to radius R 0 + na, with a a cell diameter and n ~ 2. This “forward sweep” 
has the same update rules described in section [2] for spherical range expansions. After 
the sweep, the outermost shell of the population is treated as an initial condition to 
evolve the population backwards to radius Rq — na, using the time-reversed version of 
the update rule discussed in section [2l Specifically, cells compete to divide into empty 
sites chosen in the reversed order from the forward sweep. The forward and backward 
sweeps are repeated many times, creating a dividing population of a fixed size at a 
distance P 0 from the sphere center. For more details, see |22j . Note that our results for 
how the finite size of the population influences the evolutionary dynamics at the fitness 
collapse transition should be insensitive to the particular details of how the treadmilling 
population is established. 

For M — 2, an important “rapidity reversal” symmetry exists at the DP transition 
[15] . This symmetry affects the survival probability P(t) of a mutant cluster formed 
from an initial condition with a single m — 0 class cell surrounded by cells in the rn — 1 
class. Rapidity reversal symmetry insures that, at long times, P(t) is proportional to 
the mutation-driven decay of the fraction fo(t) of 0-class cells starting from an initial 
population entirely in class m = 0. In particular, P(t ) a nf 0 {t) ~ P 5(0) at the transition, 
where h® ~ 0.451 is the DP critical exponent [T5] . The constant of proportionality 
hi depends on p and approaches zero as hi ~ yjJ1 when fi —> 0. In inflating, two- 
dimensional (circular) range expansions, this rapidity reversal symmetry is broken after 
the cross-over time t* = Rq/v (where v is the front speed) because the wandering of the 
mutant cluster boundaries gets overwhelmed by the inflating perimeter [9|. We expect a 
similar symmetry violation with respect to rapidity reversal in three-dimensional range 
expansions. 

The rapidity reversal results are shown in figure [9j We indeed find that P(t) ~ 
Kfo(t), with hi fa 0.6, at early times such that t -C t*. However, if a mutant 
cluster survives past the cross-over time t*, the cluster will inflate along with the 
frontier, thus allowing it to survive indefinitely with some limiting survival probability 
Poo = P(t —>■ oo). Our simulations indicate that Pit) saturates after time t > t *, 
approaching a constant given approximately by (t*/T g )~ s{ ° ] (black points in figure [9]). 
Since t* = Ro/v, we estimate that P Q0 ~ (t* / T g )~ &W = (a/Ro) 6ia) near the fitness collapse 
transition. Note that this result is dramatically different from a non-inflating, three- 
dimensional range expansion, for which we find P 00 = 0 at the transition. This survival 
probability enhancement is also present at the voter model point, where P^ ~ a/R q for 
spherical expansions [22). Conversely, the fraction f Q (t) of 0-class cells will continue to 
decrease at the transition when we have an initial population of m = 0 cells. So, just 
as in two-dimensional circular range expansions, rapidity reversal is broken after time 
t* in three-dimensional spherical range expansions. 

We now study other properties of the cluster formed from a single m = 0 mutant 
cell to better understand the effects of inflation. Two key quantities are the average 
squared cluster spread (X 2 (t)) and the average number of 0-class cells (. N 0 (t )) in all 
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Figure 9. Rapidity reversal symmetry in spherical range expansions spoiled by 
inflation near the fitness collapse transition (g = 0.1, s = 0.2586). The red symbols 
correspond to the decaying fraction of 0-class cells fo(t) for an initial population entirely 
in the m = 0 class. The black symbols show the survival probability P(t) of a cluster 
generated from a single 0-class cell. Both f 0 (t) and P(t) are averaged over at least 
3 x 10 4 runs. We find n « 0.6 by adjusting the data collapses until these two quantities 
overlap at early times. Note that these quantities decay in a similar way for t <C t*, 
consistent with DP scaling, indicated by the solid blue line. At late times, however, the 
survival probability approaches a non-zero constant P(t — > oo) ~ (t* /T g )~ s( °\ where 
(5(°) « 0.451 (HT, while f Q (t) continues to decay. 


surviving 0-class clusters at time t. The spread X(t) is the arc length between the 
position of the initial m — 0 cell and a rn — 0 cell in the surviving cluster at time 
t. We average X 2 {t) over all m = 0 cells at the frontier at time t and over many 
simulation runs. In the inflationary regime t > t*, the area covered by the cluster 
increases approximately quadratically in time t, due to inflation. Inflation thus leads 
to the long time scaling (X 2 (t)) ~ t 2 shown in figure fTUT a). Inside this quadratically 
increasing area, a critical directed percolation process occurs, with a decaying mutant 
fraction fo(t) ~ t~ 5< '° ) . Upon combining the quadratic area scaling with the mutant 
fraction decay scaling, we find that (. N 0 (t )) ~ ( t/t*) 2 ~ 5{0) for t t*, as shown by the 
solid black line in figure fTUT b). These results indicate that the large scale features of 
the cluster, such as its spatial spread at time t, are dictated by the inflating population 
frontier. Local features such as the local 0-class fraction decay, however, retain their 
DP properties. 

To highlight the differences, we now compare the scaling functions for inflating 
spherical expansions to those for a treadmilling population on a sphere (see figure [TTh. In 
a treadmilling population, the finite size of the population introduces strong corrections 
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Figure 10. Scaling of the cluster formed from a single mutant m = 0 cell near the 
transition. The solid lines in both panels show the expected DP scaling for t < t* and 
inflationary scaling for t t*, where z « 1.7660(16) is the dynamic exponent for DP 
m- All results are averaged over at least 3 x 10 4 runs. In (a), we show the average 
squared spread ( X 2 (t )) of all clusters that survive at least until time t. When t t*, 
inflation takes over and the mean cluster spread locks into the linear increase of the 
inflating sphere radius in time. In (b), we track the number N(t) of m = 0 cells in 
the cluster at time t. At long times 1 > i*, the cluster spread increases linearly in 
time, but the cell fraction inside the cluster decreases as f _l5( ) , yielding the long time 
scaling behavior indicated by the solid line. 


to the scaling behavior at long times. However, the corrections are quite different from 
the corrections due to inflation. For example, at long times, the one-way mutations 
lead to an even more rapid (exponential in time) extinction of the fittest m = 0 cells 
at the transition, as shown in figure □Ha). By contrast, inflation is able to save the 
m — 0 sector with a non-zero probability. Also, we see in figure fllT b) that the sector 
size saturates due to the finite population front size in a treadmilling population. In the 
inflating case, the sector grows even more rapidly at times t t*. 
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Figure 11. A comparison of (a) the survival probability and (b) the average cluster 
spread for treadmilling and inflating spherical populations. The red crosses for the 
inflating expansions in (a) and (b) are derived from the same data as in figure [9] and 
figure niil'b). respectively, with again the initial radii Rq = 5,10, 30, 50, 70, 90,110. The 
black circles are a data collapse of many treadmilling spherical populations fixed at the 
initial radii Rq of the inflating expansions. Note that the inflating and treadmilling 
expansions have different scaling variables, shown along the vertical and horizontal 
axes with the inflationary case shown first. In (a), we see that contrary to an inflating 
expansion, mutations always lead to the extinction of a m = 0 strain in a treadmilling 
population at the transition. In (b), we find that the m = 0 cluster saturates at a 
finite size related to the population front size in a treadmilling tumor. At early times, 
the treadmilling and inflating tumors have similar scaling properties. 


6. Conclusions 

Range expansions near a fitness collapse transition present a fascinating example of a 
non-equilibrium critical phenomenon. We have shown how two- and three-dimensional 
range expansions with thin, actively growing frontiers exhibit multicritical scaling 
behavior which may be characterized by coupled directed percolation processes. We 
proposed a stochastic partial differential equation hierarchy to describe these expansions 
and analyzed the hierarchy using mean-field theory. For three-dimensional expansions, 
we varied the irreversible, deleterious mutation rate /i and the strength s of the mutation 
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to calculate the shape of the phase diagram. We pointed out key differences from two- 
dimensional spatial range expansions and examined the effect of varying the possible 
maximum number of accumulated mutations, M—1 (M total species). In particular, the 
weaker effects of genetic drift yield a more stable fitness distribution in three dimensions, 
with fitness collapse occurring in a smaller region of the phase space. We also found 
that increasing M leads to possible deviations (or a slow crossover) from the expected 
directed percolation scaling predictions. 

We also considered how inflating population frontiers modify the fitness collapse 
transition for the case of just a single deleterious mutation (M = 2 total species), which 
is relevant for studying slightly deleterious passenger mutations in cancerous tissue. 
We find dramatic differences between treadmilling and inflating expansions, similar to 
those highlighted in (22] . We expect that the broad features of our results, such as the 
enhanced survival probability and cluster size scaling in the inflationary regime, will 
survive in a more realistic model with many accumulating mutations (M —>■ oo). 

There is much room for future work. For example, simulations of larger range 
expansions would be necessary to confirm that we are not seeing transient behavior and 
a universality class genuinely different from directed percolation in three-dimensional 
expansions. An exploration of the proposed equation for the evolutionary dynamics, 
equation (J3J), (via a renormalization group analysis or a similar technique) would also 
be helpful in this context. Another interesting direction would be to introduce a 
more realistic growth dynamics that allows the range expansion to develop a rough 
front. These rough fronts can couple strongly to the evolutionary dynamics and lead to 
profound changes in the behavior of the evolution near non-equilibrium phase transitions 
[251126]. 
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